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' Abstract 

The first prediction of the probability distribution function (PDF) of self-organized shear flows 
is presented in a nonlinear diffusion model where shear flows are generated by a stochastic forcing 
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coherent structure is utilized for the prediction of the strongly intermittent exponential PDF tails 
of the gradient of shear flows. Numerical simulations using Gaussian forcing not only confirm these 

"ca ■ 

>^ predictions, but also reveal the significant contribution from the PDF tails with a large population 

Qh' of super-critical gradients. The validity of the nonlinear diffusion model is then examined using a 

^ | threshold model where eddy diffusivity is given by discontinuous values, elucidating an important 

f — ' role of relative time scales of relaxation and disturbance in the determination of the PDFs. 
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I. INTRODUCTION 



Many important phenomena in nature are often far from equilibrium, strongly driven 
by instabilities or by external forces. Examples are diverse, from forest-fires to interstellar 
turbulence, which is constantly stirred by super-nova explosions. Multi-scale interactions 
are responsible for inevitably complex dynamics in these non-equilibrium systems, a proper 
understanding and description of which remains as a significant challenge in classical physics. 
As a remarkable consequence of multi-scale interactions, a quasi-equilibrium state can how- 
ever be maintained by hovering around a marginal state by a continuous adjustment of 
perturbations to establish a new equilibrium [1]. While a small perturbation can be trans- 
ported by the excitation of waves around a quiet equilibrium state, the relaxation of a large 
deviation can involve ballistic, avalanche- like events of large amplitude on a short dynamical 
time scale. This is an essential feature of the so-called self-organization, or self-organized 
criticality (SOC) in a more restricted sense [2]. In particular, it appears to be a power- 
ful paradigm for understanding complexity in plasmas, with a growing body of supporting 
evidence for self organization from computer simulations, experiments, and observations in 
laboratory and astrophysical plasmas [3, 4, 5, 6, 7, 8, 9, 10]. 

The purpose of this paper is to provide a statistical theory of self organization, which can 
perhaps be utilized as an exploratory model in different contexts. As a concrete example, we 
consider a forced shear flow, whose gradient grows until it becomes unstable according to the 
stability criterion. For instance, in a strongly stably stratified medium, fluctuations on small 
scales (or internal gravity waves) will sharpen the structure of a shear flow u [11, 12], acting 
as a forcing, until its gradient d x u = u x exceeds the critical value u xc , set by Richardson 
criterion R> R c = (u xc /Af) 2 = 1/4. Here M is the buoyancy frequency due to the restoring 
force (buoyancy) in a stably stratified medium. Once it becomes unstable, the shear flow will 
relax its gradient rapidly and generate turbulence (fluctuations) until it starts building up 
again at the expense of fluctuations. In magnetically confined plasmas, poloidal shear flows 
(zonal flows) and/or parallel flows can be generated from drift waves while becoming subject 
to Kelvin-Helmholtz type instabilities [3, 13]. Although precise physical mechanisms for the 
generation and damping may differ, the repetition of growth and damping of shear flow is 
generic, occurring in many other systems, playing a crucial role in momentum transport, 
mixing, etc. 
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We model the essential physics involved in the self-organization alluded to above by the 
following one dimensional (ID) nonlinear diffusion equation for u x [14], 



d t u x = d xx [D(u)u x ] + /, 



(1) 



where 



D(u) = v + Pu 2 x 



(2) 



In Eqs. (l)-(2), / is an external forcing; D(u) represents the effective diffusion coefficient 
including both the molecular diffusivity v and nonlinear (eddy) diffusivity capturing relax- 
ation process for unstable shear flow \u x \ > u xc . A similar quadratic eddy diffusivity has 
widely been used in modelling chemical mixing and angular momentum transport (e.g. in 
stars and the Sun) although the precise value of parameter (5 has been controversial, often 
adjusted in an attempt to reproduce observational data [15]. 



Since the nonlinear diffusion in Eq. (2) becomes important for large gradient \u x \ > J v//3, 



inhibiting further increase in the gradient, the critical gradient is roughly u xc ~ yv/P- 
Due to the relaxation of the gradient above this critical value, a naive expectation is that 
the value of gradient mostly remains subcritical. This would however be the case only 
when a relaxation time is sufficiently short compared with the characteristic time scale 
of the perturbation (forcing). In the realistic situation with continuous perturbation (a 
stochastic forcing), there will be a broad distribution of the gradient, some values of which 
may exceed far above its critical value. The key quantity to be determined is thus the 
probability distribution function (PDF) of the gradient, rather than its average value. In 
the following, we provide the first prediction of the PDF tails of the gradient by analytical 
and computational studies. Specifically, we show that PDFs tails for large \u x \ > u xc are 
strongly non-Gaussian (intermittent) with an exponential scaling exp (— cu x ) while near the 
center for small \u x \ < u XC) PDFs are Gaussian exp (—cu x ) (c = constant). We then discuss 
a threshold model where D(u) is given by discontinuous values to examine the validity of 
the nonlinear diffusion model, elucidating a crucial role of relative time scales of relaxation 
and disturbance in the PDFs. Note that there has been a growing interest in statistical 
analysis of SOC by using various statistical measures, including PDFs of avalanches [7, 8]. 

The remainder of the paper is organized as follows. We present analytical and numerical 
results of the PDFs of shear in a nonlinear diffusion model in §2. The predicted power 
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spectrum in this model is provided in §3. A threshold model is investigated in §4, with 
numerical results presented. Section 5 contains Conclusion. 



II. PDFS OF SHEAR IN A NONLINEAR DIFFUSION MODEL 

In this section, we provide analytical prediction and numerical simulation results of the 
PDFs of the shear u x in a nonlinear diffusion model (l)-(2), in particular, showing the 
agreement on strongly intermittent exponential PDF tails of u x . 

A. Analytic result 

Since for small \u x \ <C u xc , the forcing is balanced by linear diffusion, naturally leading 
to the Gaussian distribution of u x , we focus on the PDF tails for large value of \u x \ where 
the cubic nonlinearity becomes important. In order to incorporate this nonlinear interaction 
non-perturbatively, our key idea is to look for a nonlinear structure that is likely to be 
naturally sustained in a system. One candidate for such a nonlinear structure is an exact 
nonlinear solution u x oc x to Eqs. (l)-(2) in the absence of the forcing. Due to a stochastic 
forcing, this structure is then likely to form in a random fashion with the temporal behaviour 
governed by Q(t) as u x ~ iQ(t)x. Note that similar coherent structures (ramps) have also 
been successfully used in the prediction of the intermittent PDF tails of (positive) velocity 
gradient in Burgers turbulence [16, 17] that agree with numerical results. The PDF of u x 
then becomes equivalent to that of Q(t), which satisfies: 



where g is the time dependent part of the forcing with the spatial profile oc x. In the case 
of a temporally short-correlated forcing 



the Fokker-Planck equation for the PDF of Q can be derived by using a standard technique 
[19]. To this end, we introduce the generating function Z(\,t) = e~ tXQ to obtain 



d t Q = -PQ 3 + g, 



(3) 



<<7(fi)<7(f 2 )> = <f(*i-*2)G, 



(4) 



d t {Z)=[3\d xxx (Z)-\ 2 G(Z), 



(5) 
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where (Z) = J dQP{Q,t) e - iXQ = P(X,t) is Z averaged over the forcing g, which is equal to 
Fourier transform of P(Q,t). The Fourier transform of Eq. (5) then gives us the evolution 
equation for P as: 

d t P(Q,t) = POqIQ'P] + Gd QQ P. (6) 

A stationary solution of Eq. (6) can easily be found to be P(Q,t) ~ -Po ex P [~PQ 4 /(4G)], 
leading to 

P(u x ; x, t) ~ P exp [-[3u 4 x /4G] . (7) 

The PDFs tails in Eq. (7) are non-Gaussian, intermittent with the exponential scaling of P ~ 
exp (— (3u x /AG), and symmetric under the reflection x — > —x (unlike Burgers turbulence [16, 
17, 18] which is anti-symmetric). This exponential tail is one manifestation of intermittency 
caused by a coherent structure. 

To complement the Fokker-Plank approach, it is instructive to consider an alternative 
non-perturbative method based on a path integral formulation [17, 20, 21, 22]. A key 
concept in this method is similar to what was alluded in our Fokker-Plank approach in that 
a temporally localized, nontrivial vacuum state with a coherent structure - the so-called 
instanton which maximizes the path integral - causes intermittency, contributing to the 
PDF tails. Main steps involved in the computation of the PDFs by the instanton method 
are as follows . First, we express the PDFs of the velocity gradient u x — v to take the value 
of A [P(A)\ in terms of a path integral: 

P(A) = J d\e iXA - s > , (8) 

where the effective action S\ is given by 



1v 2 )v] 



S\ = —i j dxdtv[d t v — d xx (v + f3i 

\ I d x dydtv(x, t)n(x — y)v(y, t) 
+ iX J dxdtv(x)5(x - x )5(t) . (9) 

Here, v is the conjugate variable to v — u x . By using the ansatz for temporally localized 
solutions v = F{t)4> and v = /i(£)0 in Eq. (8) and then by maximizing the effective action 
S\ with respect to F and //, we obtain the equations for F and ji (with v = 0) as follows: 

8 t F - [3c 2 F 3 = -ic 3 /i , (10) 
d t n + 3/3c 2 F 2 /i = -Ac 4 <S(f) , (11) 



where 

ci = J dx(f)(x)(f)(x) , 

cic 2 = J dx^(x)d xx [(j)(x) 3 ] , 

C1C3 = J dxdy~4>(x)K(x - y)4>{y) , 

Cl c 4 = <j>(x(t = 0)) = 0(x o ) . (12) 

Since instanton v propagates forward in time and its conjugate variable v backward in time 
while the PDF is computed at t — 0, the boundary conditions on F and fi are: 

F(-oo) = 0, fi(t >0)=0. (13) 

For t < 0, Eqs. (10) and (11) give us the equation for F(t) as 

d tt F = 3c 2 (3 2 F 5 , (14) 

which can be solved with the boundary conditions F(t = 0) = F and F(t — > —00) = [Eq. 

(13)]; 

F = — (15) 

(l + 2/3c 2 F 2 t)V2- ^ i0 ^ 

To find the value of F , we integrate Eq. (11) for an infinitesimal time interval t = [— e, e] 
by using Eq. (13) as: 

/i(-e) = /i(0) = Ac 4 , (16) 
and substitute Eq. (16) and d t F = — PC2F 3 [from Eq. (15)] in Eq. (10) to obtain 

* = W S « A ' (17) 

where q = ic 3 C4/2c2/3. 

The determination of the PDFs of P(A) now requires a few more steps. First, we evaluate 
S\ in Eq. (8) by using Eqs. (15), (16) and (17): 

Sx = QX 4/3 , (18) 

where Q = 3i(ciC 4 )g 1 / 3 /4. The next step is to find the PDF tails by computing the A integral 
(8) in the limit of large A. To this end, we substitute Eq. (18) into Eq. (8) and approximately 
evaluate A integral as / d\e iXA ~ s * = J d\e~ G{x) ~ e - G ^°\ where G(X) = -iXA + S x and 



A = (3iA/AQ) 3 is a saddle-point which minimizes G(X). Therefore, P(A) in Eq. (8) 
becomes 



(c 2 ci)ci 



(19) 



(c 3 ci) 

Equation (19) is the exponential tail of PDF of u x to take the value of A, with the same 
exponential scaling as in (7). Of importance to notice is that the result (19) follows from 
the order of the highest cubic nonlinearity in Eq. (1), being independent of the precise form 
of the spatial structure of and 0, which has not been specified yet. The latter however 
plays a crucial role in the determination of the overall amplitude of the PDFs through 
the coefficient £ [see Eq. (19)]. Fortunately, the form of oc x and oc x 3 (i.e. the exact 
nonlinear solutions to v and v) can be inferred from the instanton equations d t v — d xx {(5v 3 ) = 
—i J dyK,(x — y)v(y,t) and d t v + d xx (3(3v 2 v) = —Xv(xo)5(t), obtained by minimizing S\ with 
respect to v and v, and by then using k(x — y) ~ k q [1 — (x — y) 2 /2 + ■ • •] = k [1 + xy + • ■ ■]. 
The use of oc x and oc x 3 in Eqs. (12) and (19) then gives |c 2 Ci/c 3 | ~ 6//t , and thus 
£~3(/3/ko). 

To summarize, both Fokker-Planck and instanton methods, based on the key idea that 
the PDFs tails are caused by a coherent structure, give us the strongly intermittent PDFs 
tails of shear gradient u x , with the same exponential scalings exp (—cu\) (c = constant) [see 
Eqs. (7) and (19)]. It is interesting to compare these results with the right PDFs of the 
velocity gradient in Burgers turbulence, which was predicted to be exponential with a dif- 
ferent exponent [i.e., exp (—cu 3 )} due to ramp-like coherent structures (u oc x) [17] (followed 
by numerical verification). This scaling with the different exponent basically results from 
the quadratic highest nonlinear interaction in Burgers turbulence, different from the cubic 
highest interaction in our model (l)-(2) (see [21] for more details). In plasma turbulence, 
exponential PDFs tails of various fluxes have been theoretically predicted without numerical 
confirmation (e.g. see [20, 21, 22]). Nevertheless, it is very interesting that these exponential 
scalings have often been observed in the tails of fluxes in laboratory plasmas (e.g. see Refs. 
[23, 24]). 
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B. Numerical Results 



To test our analytical prediction (7) and (19), we perform direct numerical simulations by 
numerically integrating Eqs. (1) and (2) using method outlined in [14]. To briefly recap, we 
use finite difference method to solve (1). The spatial discretization is second order accurate 
and the time integration uses Euler-Maruyama method. Adaptive time stepping is also 
used for numerical stability of the diffusion term. For each step of the simulations, the 
Gaussian noise is produced using the Box-Muller method [25], which gives homogeneous, 
and temporally short-correlated forcing / in Eq. (1) with the power spectrum F{k): 

(f(h, h)f(k 2 , h)) = 8(h - t 2 )8(h + k 2 )F(k) . (20) 

The results for the PDFs of u x , P(u x ), for a white-noise F(k) = k° are shown by the solid 
line in Fig. 1 for the values of parameters v = 6 x 10~ 3 and /3 = 6.25 x 10~ 3 . It can clearly be 
seen that the PDF is Gaussian near the center but becomes exponential exp (—cu x ) in the 
tails (c = constant). These exponential tails agree perfectly with our theoretical prediction 
(7). To highlight this, the dotted and dashed lines in Fig. 1 are fits to a Gaussian and to 
exp (— cu x ), respectively. The cross-over between these two regimes occurs approximately at 
the expected critical gradient of u xc ~ \Jv/P = 0.98. The mean value of \u x \ is found to be 
smaller than this, with the value about 0.59. However, there is yet a significant probability of 
20% of super-critical gradient \u x \ > \u xc \ from the PDF tails. The intermittent occurrence 
of super-critical gradients can be appreciated from the profile of u x plotted in Fig. 2, which 
exhibits a bursty of large gradients. Reflectional symmetry of the PDF is also seen in Figs. 1 
and 2. 

While mathematically, the PDF tails exp (— cu x ) result from the highest cubic nonlinearity 
in the equation for u x (1), physically, they are due to the feedback of shear on turbulence 
when it becomes unstable. That is, while shear is generated by turbulence [modelled by the 
forcing / in Eq. (1)], it feeds back on turbulence, limiting its own growth, thereby reducing 
the PDF tails below the Gaussian prediction (see Fig. 1). We have confirmed that these 
exponential PDFs tails are robust features by using different power spectra F(k) ~ k~ l and 
k~\ 
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III. POWER SPECTRUM IN A NONLINEAR DIFFUSION MODEL 



One of the main interests in the previous studies of self-organization (or SOC) has been 
power spectrum. It is thus interesting to examine what prediction can be made on the power 
spectra in our model. To this end, we compute the PDFs of u x (k), i.e. P(k,t) by observing 
that the evolution of k mode involves the cubic nonlinearity due to nonlinear diffusion in 
Eq. (2) of the form / dk\dk 2 u x {ki)u x {k2)u x {k — k± — k 2 ). We approximate the latter as 
\u x {k)\ 2 u x {k) by keeping only the dominant coherent interaction (which can be justified for 
a narrow spectrum), and rewrite Eq. (1) as follows: 

d t u x ~ j3k 2 u x - vk 2 u x + / . (21) 

The Fokker-Planck equation for the P(u x ; k,t) can be obtained as previously, from which a 
stationary PDF follows as: 

P(U X - k, t) = p ( k ) e -k 2+r >4(k 2 ul/2+„) ^ (22) 

where Po(k) = 1/ Jdu x P(u x ;t) is the normalization constant, and the power spectrum 
F(k) = k~ r is used. In the linear case, it is easy to see that the power spectrum p(k) = 
(\u x (k)\ 2 ) = j P(u x ; t)\u x (k)\ 2 oc k~ r ~ 2 . On the other hand, in a strongly nonlinear case, we 
find that 

p(k) oc k~( 2+r ^ . (23) 

Remarkably, the prediction (23) agrees very well with the numerical results shown in Liu [14] . 
In particular, in the case of the red noise with r = 2, Eq. (23) predicts k~ 3 spectrum, with a 
better agreement with numerical result than the prediction (/c~ 3 5 ) from the renormalization 
theory! 

IV. THRESHOLD MODEL 

Our results highlight the importance of the statistical description of self-organization. In 
particular, the population of super-critical gradients as well as the form of PDF tails can 
depend on the relative time scales between disturbance (i.e. forcing) and relaxation. To 
show this, we consider a threshold model where the nonlinear diffusion D{u) in Eq. (1) is 



9 



given by the two discrete values as 

[v, for \u x \ < u xc ; 

D{u) = 

[V (> v), for \u x \ > u xc . 

Here v is molecular diffusivity while V represents a large diffusion due to avalanche- like 
events which efficiently relax super-critical gradients [14]. To investigate the extreme limit 
where the relaxation rapidly occurs on the shortest time scale, we numerically solve Eq. (1) 
by using this discrete D{u) and by applying forcing when \u x \ is less than u xc everywhere in 
the domain in order to ensure that the relaxation occurs much faster than the disturbance. 
Note that a similar method was used, for example, in the SOC solar flare models [9, 10]. 
Numerical simulation results using v = 6 x 10~ 3 , V = 4.5 x 10~ 2 , and 

W^c — 2 £1X6 plotted 

in Fig. 3, which shows that only 0.24% shear is super-critical. This is much less than 
20% found in Fig. 1 in the nonlinear diffusion model, and is due to rapid relaxation by a 
large diffusion V for \u x \/u xc > 1. The resulting PDFs for this gradient \u x \/u xc > 1 are 
Gaussian as can be seen in Fig. 3 since the diffusion in Eq. (1) is essentially linear. In 
comparison, the Gaussian PDFs near the center for small \u x \/u xc < 0.34 results from small 
fluctuations which satisfy Gaussian statistics. What is very interesting is that there is a 
window of piece-wise exponentials exp (—cu x ) between these two Gaussian PDFs for the 
gradient 0.34 < \u x \/u xc < 1, with a significant population 30%. This exponential PDFs 
are similar to those found in the nonlinear diffusion case for \u x \/u xc > 1 although the 
exact values of \u x \/u xc for the exponential PDFs are not identical. Therefore, these results 
indicate that the nonlinear diffusion can be a reasonable approximation for a certain range 
of the shear values and relaxation time scales. 

We have also performed the simulation by applying both the forcing and diffusive relax- 
ation simultaneously to make disturbance time sufficiently short. The resulting PDFs are 
found to be Gaussian since the diffusion [Eq. (1)] in this case is essentially linear except at 
\u x \ = u xc . The super-critical population is also higher (19%) due to the slower relaxation. 



V. CONCLUSION 



We have presented a statistical theory of self-organization by utilizing a simplified nonlin- 
ear diffusion model for a shear flow and a widely invoked quadratic eddy diffusivity [14, 15]. 
Both Fokker-Planck and instanton methods predict the PDF tails of the exponential form 
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exp (— ci4), with a strong intermittency. Our numerical simulation using Gaussian forcing 
with three different power spectra not only confirm these predictions, but also reveal the sig- 
nificant contribution from the PDF tails with a large population of super-critical gradients, 
which could play a crucial role. These results highlight the importance of the statistical 
description of gradients in self-organization, rather than its average value as has conven- 
tionally been done. The validity of the nonlinear diffusion model was then examined using 
a threshold model, elucidating an important role of relative time scales of relaxation and 
disturbance, calling for a care in actual modelling of a particular system. 

Our results can have significant implications for the dynamics and the role of shear 
flows (e.g. zonal flows) in laboratory, astrophysical and geophysical plasmas, which is vital 
not only in momentum transport, but also in transporting chemical species and controling 
mixing of other quantities (e.g. air pollution, weather control) [26, 27]. Our theory 
can also provide a useful guide in understanding self-organization in other disciplines, 
such as population in environmental dynamics and biology, forest-fire, and reaction and 
diffusion in chemistry. Future work will include specific applications to those systems, 
the extension of our model to incorporate the finite correlation time of the forcing and 
a non-diffusive flux [28], and the investigation of the joint PDFs of fluctuations and 
mean gradients in a consistent way. Note that an initial attempt to the prediction 
of the joint PDFs has been made in the ion temperature gradient turbulence (for mag- 
netically confined plasmas) by neglecting the feedback of shear flows on fluctuations [22, 29]. 

This research was supported by the EPSRC grant EP/D064317/1 and RAS Travel Grant. 
The National Center for Atmospheric Research is sponsored by the NSF. 
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Figure Captions 



Fig. 1 The solid line is the PDFs from the numerical simulation of a nonlinear diffusion 
model (2) for a white noise. The dotted and dashed lines are the fits to Gaussian and 
exp (— ci4) (c =const). 

Fig2. The profile of u x corresponding to Fig. 1. 

Fig. 3 Solid line is the PDF from a threshold model. Dotted and dash-dotted-dotted-dotted 
lines are Gaussian fits; dashed and dashed-dotted lines are fits to exp(— cu^) (c =const). 
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